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Abstract. We present the results of large scale computer simulations in which 
we investigate the structural and dynamic properties of silicate melts with the 
compositions (Na20)2(Si02) and (Al203)2(Si02). In order to treat such systems 
on a time scale of several nanoseconds and for system sizes of several thousand 
atoms it is necessary to use parallel supercomputers like the CRAY T3E. We show 
that the silicates under consideration exhibit additional intermediate range order 
as compared to silica (Si02) where the characteristic intermediate length scales 
stem from the tetrahedral network structure. For the sodium silicate system it is 
demonstrated that the latter structural features are intimately connected with a 
surprising dynamics in which the one-particle motion of the sodium ions appears 
on a much smaller time scale than the correlations between different sodium ions. 



1 Introduction 

Silicate melts and glasses are an important class of materials in very different 
fields, e.g. in geosciences (since silicates are geologically the most relevant 
materials) and in technology (windows, containers, optical fibers etc.). From 
a physical point of view it is a very challenging task to understand the prop- 
erties of those materials on a microscopic level, and in the last twenty years 
many studies on different systems have shown that molecular dynamics com- 
puter simulations are a very appropriate tool for this purpose The 
main advantage of such simulations is that they give access to the whole 
microscopic information in form of the particle trajectories which of course 
cannot be determined in real experiments. 

In pure silica (Si02) the structure is that of a disordered tetrahedral net- 
work in which SiC>4 tetrahedra are connected via the oxygens at the corners. 
In recent simulations we have studied in detail various aspects of static and 
dynamic properties of silica such as structural and thermodynamic prop- 
erties of theglass state |||(|, the diffusion dynamics and structural relax- 
ation fTjjBpyifjf , the frequency dependent specific heat the vibrational 
degrees of freedom ]l2| and free surfaces Jl3|,[l4 15 . In this paper we con- 



sider silicates that contain additional oxide components. Especially silicates 
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with alkali oxides have been investigated very recently in several molecular 
dynamics simulations [l^ . p^| , p^| , ^9pc|pl| , p^ | . We investigate here the systems 
(Na 2 0)2(Si0 2 ) and (Al 2 3 )2(Si0 2 ) (for which we use in the following the 
abbreviations NS2 and AS2, respectively). Whereas sodium in NS2 plays 
the role of a network modifier that partially destroys the SiC>4 network, alu- 
minium in AS2 is also a network former in that it prefers a four-fold co- 
ordination by oxygen atoms. However, the packing of the AIO4 tetrahedra 
is different from that of the SiC>4 tetrahedra which is indicated for instance 
by a different coordination number distribution of aluminium and silicon 
by oxygen atoms (mainly two-fold for silicon and two- and three-fold for 
aluminium) . As we show in the following the insertion of sodium or alu- 
minium atoms does not only modify the structure on local length scales but 
it introduces also new intermediate length scales that can be identified by 
means of the partial static structure factors. These length scales are impor- 
tant for the dynamic properties as we will demonstrate for the case of NS2. 

The rest of the paper is organized as follows: In the next section we give 
the main computational details and discuss the efficiency of our simulation 
code on the T3E at the HLRZ Stuttgart. Then we present the structural 
properties of AS2 and NS2 on intermediate length scales (Sec. 3) and the 
dynamics of NS2 (Sec. 4). Eventually we summarize our results (Sec. 5). 

2 Model and details of the simulations 

In a classical molecular dynamics (MD) computer simulation one solves nu- 
merically Newton's equations of motion for a many particle system. If quan- 
tum mechanical effects can be neglected such simulations are able to give 
in principle a realistic description of any molecular system. The determining 
factor of how well the properties of a real material are reproduced by a MD 
simulation is given by the potential with which the interaction between the 
atoms is described. The model potential we use to compute the interaction 
between the ions in NS2 and AS2 is the one proposed by Kramer et al. 
which is a generalization of the so-called BKS potential ]25| for pure silica. 
It has the following functional form: 

^ r) = ?^l + A af3 exp(-B a0 r)-^f a,(3 £ [Si, Al, Na, O]. (1) 

Here r is the distance between an ion of type a and an ion of type (3. The 
values of the parameters A a p, B a p and Cap can be found in the original pub- 
lication. The potential (|l|) has been optimized by Kramer et al. for zeolites, 
i.e. for systems that consist of Si, Al, Na, O, and possible other components 
like phosphor. In that paper the authors used for silicon, aluminium, and 
oxygen the partial charges gsi = 2.4, <jai = 1-9, and qo = —1.2, respectively, 
whereas sodium was assigned its real ion charge <?Na = 1-0. Thus, with this 
set of charges charge neutrality is fulfilled neither in NS2 nor in AS2. We 
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number of processors 

Fig. 1. Speed up factor for the simulations with N — 8064 (filled triangles) 
and N = 1408 particles (open circles) as a function of the number of pro- 
cessors n. The bisecting line (straight line) indicates a perfect scaling of the 
performance with n. 



have therefore modified the Kramer potential by setting the partial charge 
for sodium and aluminium to 0.6 and 1.8, respectively, and by introducing ad- 
ditional short range potentials such that the original functional form of the 
Kramer potential is approximately recovered on distances of nearest Al-0 
and Na-0 neighbors. More details on the interaction potential can be found 
in Refs. Jl^,[l7|,^3) . Our models give predictions for structural and dynamic 
properties of NS2 and AS2 which are in good agreement with experimental 
findings |l7| , ^3| . Furthermore, Ispas et al. [^6| have shown for (Na20)4(SiC>2) 
that ab initio simulations (Car-Parrincllo molecular dynamics) yield compa- 
rable results regarding the structure to those obtained with molecular dy- 
namics simulations in which our potential model is used. 

The simulations have been done at constant volume: For AS2 we fixed the 
density to 2.6 g/cm 3 which is close to the experimental density at T — 300 K. 
In the case of NS2 we did simulations at the two densities 2.37 g/cm 3 and 
2.5 g/cm 3 , corresponding to experimental densities in the melt and at room 
temperatures, respectively. The AS2 system consists of 1480 particles and for 
the NS2 systems we used system sizes of 8064 particles at the low density 
and 1152 particles at the high density. 
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As can be seen from Eq. [l] the interaction potential contains a long-ranged 
Coulomb term. This part of the interaction is the most time consuming in 
the calculation of the forces. To do this we made use of the so-called Ewald 
summation technique p7| , a method that scales with the particle number 
N as TV 3 / 2 . Thus, for systems which contain about 8000 particles a huge 
numerical effort is required: The longest runs (at the lowest temperatures) 
had a length of about 10 million time steps for which a time of two weeks was 
needed on 64 processors thus giving a total CPU time of about 128 weeks of 
(single) processor time. 

The equations of motion were integrated with the velocity form of the 
Verlet algorithm. The time step of the integration was 1.6 fs. The temper- 
ature range investigated was 4000 K> T > 2100 K in the case of NS2 and 
6100 K> T > 2300 K in the case of AS2. To equilibrate the systems the 
temperatures were controlled by coupling them to a stochastic heat bath, i.e. 
by substituting periodically the velocities of the particles with the ones from 
a Maxwell-Boltzmann distribution with the correct temperature. After the 
system was equilibrated at the target temperature, we continued the run in 
the microcanonical ensemble, i.e. the heat bath was switched off. We have 
done production runs up to several ns real time which corresponds to several 
million time steps. We have also calculated glass structures at T = 300 K. 
The glass state was produced by cooling the system from equilibrated con- 
figurations at our lowest temperatures with a cooling rate of about 10 12 K/s. 
Note that we show in the following sections only results for the lowest tem- 
peratures, i.e. at T = 2100 K for NS2 and at T = 2300 K for AS2 as well as at 
T = 300 K for both systems, because the results for the higher temperatures 
lead essentially to the same conclusions that we will draw below. However, 
a detailed discussion of the temperature dependence of the systems under 
consideration can be found in Refs. [|17j . 

The program code was written in FORTRAN. All the parallelization was 
done by using MPI subroutines. More details on the parallelization can be 
found in Refs. pp3]|. Of course, the performance of a parallel code never 
scales perfectly with the number of processors because the communication 
between the processors requires an additional amount of CPU time. Fig. [l] 
shows the speed up factor on the Cray T3E of the HLRZ Stuttgart as a 
function of the number of processors n, i.e., the factor by which the code 
is faster if one uses n processors instead of one. The bisecting line indicates 
the limiting case where the communication overhead is not influenced by the 
speed of the code. We see that the curves for N — 1408 and N — 8064 scale 
nearly perfectly for n < 16. For n — 64 we obtain still a speed up factor of 
about 46.4 for N = 8064 particles whereas this factor is 42 for N = 1408. In 
most of our simulations we have used 64 processors for the large systems and 
32 processors for the small ones. 
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Fig. 2. Partial static structure factors for AS2 at T = 2300 K. a) Al-Al, 
Si-Si, and Si-Al correlations, b) Si-O, Al-O, and 0-0 correlations. For the 
meaning of the dashed vertical lines see text. 



3 Intermediate length scales in silicates 

An appropriate quantity to investigate the structure of atomic systems on 
intermediate length scales is the static structure factor. It is essentially the 
Fourier transform of the pair correlation function which gives the probability 
of finding an atom at a distance r from another atom The structure 
factor can be directly measured in neutron scattering experiments from the 
intensity of the radiation observed with a momentum transfer fi,q (h: Planck's 
constant, q: wave-vector of the momentum transfer). In a three-component 
system one can define six partial structure factors as Q 

s Q/3 (?) = ^ E E < ex p fa ■ [>■* - r *])> ■ ( 2 ) 

k=l 1=1 

where the first sum runs over N a particles of type a and the second one over 
Np particles of type (3. 

Fig. | shows S a p(q) for AS2 at the temperature T = 2300 K. For q > 
2.3 A -1 the partial structure factors reflect length scales of nearest neigh- 
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Fig. 3. Snapshot of (Al 2 3 )2(Si0 2 ) (AS2) at T = 300 K. The size of the 
spheres is chosen such that one can identify aluminium- and silicon-rich 
regions: The aluminium and oxygen atoms are shown respectively as big blue 
and gold spheres, whereas the oxygen atoms are shown as small red spheres. 



bors. In AS2 the smallest distances between atoms are those of Al-0 and 
Si-0 bonds that have lengths of about 1.6 to 1.65 A. The peaks around 
</2 = 1.7 A -1 in S a p(q) (marked by dashed vertical lines in Fig. ||) are due 
to the order that arises from the tetrahedral network structure. The length 
scale 27r/1.7 A -1 = 3.7 A that corresponds to this peak is approximately 
the spatial extent of connected SiC>4 and AIO4 tetrahedra. Note that in sil- 
ica a peak at 1.7 A" 1 is also a very prominent feature and is called there 
first sharp diffraction peak. But in contrast to silica one observes in AS2 an 
additional peak at qi ~ 0.5 A -1 in the Al-Al, Si-Si, and Si-Al correlations 
and only weakly pronounced also in the remaining correlations in which oxy- 
gen is involved. q± corresponds to a length scale of about 12.5 A and has its 
reason in a slightly different ordering of AIO4 complexes as compared to the 
SiC>4 network (for details see [^3]). This relatively large length scale shows 
that large system sizes are required to analyze the structure of systems like 
AS2 in a sensible way. The different ordering of AIO4 leads to a structure 
where an AIO4 tetrahedron prefers to be surrounded on a local scale by other 
AIO4 tetrahedra. This leads to a structure where AIO4 complexes are con- 
nected to each other as string like objects through the system that form a 
percolating network. This is illustrated by the snapshhot in Fig. || where the 
aluminium and silicon atoms are shown as the blue and gold spheres, respec- 
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Fig. 4. Partial static structure factors for NS2 at T — 2100 K. a) Na-Na, 
Si-Na, and Na-0 correlations, b) Si-O, Si-Si, and 0-0 correlations. For the 
meaning of the dashed vertical lines see text. 
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tively. Note that it does not matter that this snapshot is at T = 300 K and 
not at T = 2300 K as the structure factors in Fig. |^ because we find only 
small differences in structural quantities at both temperatures. Thus, we see 
that the aluminium atoms are not at all homogeneously distributed and if 
one only considers the Al atoms voids with a size of about 2n/qi are found 
that lead to the peak at q\ in Sai-ai- It is not surprising that these voids are 
also reflected in the Si-Si and Si-Al correlations but much less in the corre- 
lations containing oxygen (as can be seen in Fig. ||b): The oxygen atoms are 
essentially homogeneously distributed on the length scale 2n/qi since they 
are nearest neighbors of silicon and aluminium with a similar length of Si-O 
and Al-0 bonds. 

The sodium ions in NS2 play a different role from the aluminium atoms 
in AS2 since they partially destroy the Si04 network. This can be directly 
recognized in the partial structure factors for NS2 which are shown in Fig. |^ 
at T = 2100 K for the density p = 2.37 g/cm 3 : The peak at q 2 = 1.7 A" 1 that 
reflects the structure of a tetrahedral network is absent in the correlations 
with sodium (Fig. ^a) and is especially in Ssi-si much weaker pronounced 



8 J. Horbach et al. 




0.0 1.0 2.0 3.0 4.0 5.0 6.0 
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Fig. 5. S n (q) at T = 300 K for AS2 and for NS2 at the indicated densities. 
The dashed vertical lines mark the position of the peaks at q\ = 0.5 A" 1 in 
AS2, qi = 0.95 A -1 in NS2, and q 2 = 1.7 A" 1 in both systems. 



than in AS2 (Fig. [|b). But we observe again a second prepeak at smaller q, 
now around q\ — 0.95 A -1 . This q value is of the order of the length scale of 
next nearest Na-Na or Si-Na neighbors (around 6.6 A). Again, the peak at q\ 
is the characteristic wave-vector of a percolating network that is now formed 
by the sodium atoms. At first glance it seems to be surprising that also Sb-O 
exhibits a peak at q%. But the role of oxygens is different in NS2 from that in 
AS2: The nearest neighbor distance for Na-O, 2.2 A, is larger than for Si-0 
which is 1.6 A. And the arrangement of oxygen around sodium is different 
from the tetrahedral one around silicon (for more details see Ref. [|l7|). Thus, 
the distribution of oxygen atoms in NS2 is not homogeneous on the length 
scale 2n/qi. 

So far we have seen that NS2 and AS2 exhibit intermediate order on a 
relatively large length scales. This gives rise to a prepeak in S a/ 3(q) at qi 
which is 0.5 A" 1 for AS2 and 0.95 A" 1 for NS2. But does one see these 
peaks at q\ also in experiments? In experiments such as neutron scattering 
one does not have access to the partial structure factors for systems like NS2 
or AS2. Here one measures a sum of the partial structure factors whereby the 
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different contributions are weighted by the neutron scattering lengths b a : 

Sn ^ = V- \ T h 2 bkbl ^ XP (* q ' ^ ~ r 'l^ ' 

The values for b a are 0.4149-10" 12 cm, 0.3449-10~ 12 cm, 0.363-1CT 12 cm, and 
0.5803-10~ 12 cm for silicon, aluminium, sodium, and oxygen, respectively p8| . 
By weighting the S a p (q) from our simulation with the b a in accordance with 
Eq. (||) one can easily calculate the quantity S n (q). It is shown in Fig. || at 
T = 300 K for NS2 at the two densities p = 2.37 g/cm 3 and 2.5 g/cm 3 and 
for AS2 at p = 2.6 g/cm 3 . We infer from this figure that the aforementioned 
prepeaks at q\ can be seen in AS2 and in NS2 at the higher density only as 
a weakly pronounced shoulder. Thus it would be difficult to identify them in 
a neutron scattering experiment. Only in NS2 at the lower density one can 
clearly see the prepeak at q = 0.95 A -1 . But at this density we observe a 
negative pressure of about —1.6 GPa at T = 300 K, a condition that would be 
difficult to realise in an experiment. However, in an experiment under normal 
pressure conditions the density decreases if one goes to higher temperatures. 
And indeed, very recent neutron scattering experiments of Meyer et al. do find 
the feature at q\ for NS2 (29). Meyer et al. have measured for the first time 
the temperature dependence of the structure factor from T = 300 K (where 
the system is in a glass state) to T — 1600 K (where one has a melt). They 
find that the feature at q\ becomes more and more pronounced by increasing 
the temperature and one can clearly identify it at T = 1600 K. This behavior 
is similar to what we see in our simulations and can be understood by an 
decreasing density in the experiment if the temperature is increased. 

4 The dynamics of NS2 

In a recent simulation we have demonstrated that the dynamics in in NS2 
is much faster than the one in pure silica |L7|. Even at a relatively high 
temperature of T — 2750 K the diffusion constants of silicon and oxygen are 
two orders of magnitude larger in NS2 than in SiC>2. Furthermore, in NS2 
the sodium diffusion decouples more and more from the silicon and oxygen 
diffusion such that at temperatures T < 2500 K the dynamics of the Na 
atoms is about two orders of magnitude faster than the one of the oxygen 
and silicon atoms Jl7| . This is in qualitative agreement with the expermental 
fact that NS2 is an ion conducting material. 

Thus, since essentially the Si and O atoms do not move with respect 
to the movement of the Na atoms one may expect that sodium diffusion is 
restricted to a small subspace of the configuration space. The Si and O atoms 
form a quasi-frozen matrix for the Na atoms and it would be surprizing if 
the sodium atoms are able to diffuse into this matrix. In order to check this 
idea we have calculated a (coarse grained) probability of finding no sodium 
atom at a given location in space. Following the approach of Jund et al. [[l9| 
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Fig. 6. Swiss cheese structure factor S sc (q, t) at the indicated times. The inset 
shows the probability P(t) (see text). The circles on the curve for P(t) are 
at the times at which S sc (q,t) is shown. 



we calculate this probability by dividing the system into 48 3 cubes (of length 
L/48 ps 1.01 A). Then we calculate the probability P(t) that a cube which 
does not contain a sodium ion at time zero is also not visited by a sodium 
ion until a later time t. The time dependence of P(t) is shown in the inset of 
Fig. |^. From this graph we recognize that after 2.5 ns, i.e. after more than the 
a— relaxation time of the matrix p7| , more than 50% of the cubes have not 
yet been visited by a sodium atom. (We mention that after this time the mean 
squared displacement of the Na atoms is more than (45 A) 2 , which shows that 
these atoms have moved a large distance. On this time scale also the local 
structure of the Si-0 matrix is partially reconstructed [JlTj.) Hence we can 
conclude that on this time scale the sodium free region forms a percolating 
structure around a network of channels, i.e. it has somewhat the structure of 
a Swiss cheese. In order to investigate the structure of this percolating region 
we define a "Swiss cheese" structure factor S sc (q,t) as follows: We assign to 
each cube which has not been visited by a sodium atom until time t a point 
and we compute the static structure factor of N sc (t) — P(i)(48 3 — N^ a ) 
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t[ps] 

Fig. 7. Coherent intermediate scattering functions F(q, t) for sodium-sodium 
correlations (bold solid lines) and incoherent intermediate scattering func- 
tions F s (q,t) (dashed lines) at T — 2100 K for the indicated values of q. 



points: 

1 N BC (t) 

Ssc(q,t) = — — ( ex P(^ ' ( r k - *■!))> • (4) 

sc ^ ' k,l=l 

This quantity is shown in Fig. |^ for four different times: t = 0.56 ps, 7.7 ps, 
164 ps, and 2.13 ns. We see that S sc (q,t) has a pronounced peak at q% = 
0.95 A -1 which is also a prominent feature in <S , Na-Na('z), as we have seen 
in the preceding section. Hence we can now conclude that the peak at qi in 
^Na— NaC?) corresponds to the typical distance between the channels. Note 
that with increasing time the height of this peak increases quickly. However, 
it is clear that the peak at q\ decreases again on the time scale on which 
the matrix starts to reconstruct itself significantly and thus rearranges the 
channel structure. 

We address now the question how the sodium ions relax inside the chan- 
nels. An appropriate quantity to investigate this point are time dependent 
density-density correlation functions, i.e. the coherent intermediate scatter- 
ing function F{q, t) and its self part, the incoherent intermediate scattering 
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function F s (q,t) 0. In Fig. fj] we show F(q,t) for the Na-Na correlations 
(solid lines) as well as F s (q,t) for the sodium atoms (dashed lines) for three 
different wave-vectors: q = 0.94 A -1 , 2.0 A -1 , and 2.75 A -1 . From this fig- 
ure we infer immediately a surprising result: At q = 0.94 A -1 , i.e. at the 
characteristic q value of the sodium channel structure, F(q,t) decays on a 
time scale which is two orders of magnitude larger than the one for F s (q, t) . 
Such a strong difference cannot be explained by a de Gennes narrowing argu- 
ment Q . Instead this result can be rationalized by the idea that the sodium 
atoms move quickly between preferential sites, since this type of motion gives 
rise to a fast decorrelation of the incoherent function whereas it does not af- 
fect the coherent one. Only on the time scale of the relaxation of the Si02 
matrix also the coherent function starts to decay. Note that the slow decay 
of F(q,t) is found only for wave-vectors around 0.95 A -1 . For different q 
the function decays significantly faster as can be seen from the other curves 
shown in Fig. [?]. 

More details on the issues discussed in this section can be found in 
Ref. [L§. 

5 Summary 

Large scale molecular dynamics computer simulations as the ones presented 
in this paper for AS2 and NS2 require the use of parallel computers such as 
the Cray T3E. Only then it is possible to simulate these systems on a scale of 
several ns for system sizes which are big enough to study the structure and 
dynamics on intermediate length scales (up to 8000 particles in our case). 
Although no neutron scattering experiments can be done yet for temperatures 
above 2000 K, the structural and dynamical properties are already present 
at these high temperatures and thus, one can gain insight into features that 
one observes in experiments. Furthermore, this insight is much more detailed 
in a MD simulation than in an experiment since one has access to the full 
microscopic information in form of the particle trajectories. 

We have exploited this fact for the case of AS2 and NS2 by showing 
that these systems exhibit intermediate range order on length scales that 
are larger than the one given from the tetrahedral network structure in pure 
silica. The reason for this is a different ordering of Al-0 and Na-0 complexes 
and leads to a percolating network of these structural elements through the 
Si04 network. We have shown for the example of NS2 that this intermediate 
range order is also important to understand the dynamics: In NS2 the sodium 
ions that move through channels in the Si-0 matrix and the structure of 
these channels is connected with the prepeak in the static structure factor at 
0.95 A -1 . The presence of these channels leads to a surprising decoupling of 
the fast (single particle) sodium motion from correlations between different 
sodium atoms that decay on the time scale of the channel relaxation. 
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